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Ionic solids and melts are compounds in which the interactions are dominated by electrostatic 
effects. However, the polarization of the ions also plays an important role in many respects as 
has been clarified in recent years thanks to the development of realistic polarizable interaction 
potentials. After detailing these models, we illustrate the importance of polarization effects on 
a series of examples concerning the structural properties, such as the stabilization of particular 
crystal structures or the formation of highly-coordinated multivalent ions in the melts, as well as 
the dynamic properties such as the diffusion of ionic species. The effects on the structure of molten 
salts interfaces (with vacuum and electrified metal) is also described. Although most of the results 
described here concern inorganic compounds (molten fluorides and chlorides, ionic oxides...), the 
particular case of the room-temperature ionic liquids, a special class of molten salts in which at 
least one species is organic, will also be briefly discussed to indicate how the ideas gained from the 
study of “simple” molten salts are being transferred to these more complex systems. 


INTRODUCTION 

Molecular dynamics now is an indispensable tool for 
solving many physical chemistry problems. Alongside 
the constantly increasing efficiency of computers, the de¬ 
velopment of new modeling strategies as well as efficient 
algorithms has opened the way for a new approach to the 
study of many systems. Simulation may be employed as 
a predictive tool that can be used to complement, or 
even anticipate experimental measurements. Ideally, to 
be completely predictive, atomic-scale simulations should 
be based entirely on parameter-free quantum chemical 
calculations. This is not yet achievable, so that con¬ 
structing better models for interatomic interactions re¬ 
mains an important goal, especially for condensed mat¬ 
ter systems. |Tj That is, to reach the time scales and 
system sizes which are necessary to understand complex 
phenomena such as transport mechanisms in the liquids 
or biological processes, the use of explicit interaction po¬ 
tentials derived from a model is necessary. [5] Computer 
simulations are a particularly important predictive tool 
for molten salts, for which experiments are often difficult 
or even impossible because of the extreme physical condi¬ 
tions and because many melts are highly corrosive. EH] 
The objective for model development is a transferable 
model, i.e. one in which the corresponding interaction 
potential parameters have to be determined once and for 
all before to be used in a variety of thermodynamic con¬ 
ditions (temperature, pressure, composition). 

Halide melts have been the testing ground for this 
model development, as they were the target of pioneering 
experimental studies to examine the atomic scale struc¬ 
ture by spectroscopy EH3 and diffraction [SHU! ; more 
recently attention has shifted to oxides for which the 
experimental difficulties are more pronounced. Due to 
the monatomic nature of the halide melts, the complex¬ 
ities of intramolecular interactions and molecular shape 
are avoided both in constructing the interaction model 


and in extricating the influence of intermolecular ef¬ 
fects from experimental data. Although, as we will see, 
complex local coordination structures are present in the 
melts, mm they arise from the interplay of the atomic 
interactions. Most of the short-range structural proper¬ 
ties can be attributed to the competition between the 
overlap-repulsion on the one hand and the Coulombic 
interaction on the other hand m- The former, which 
is a consequence of the Pauli principle, will determine 
the closest approach distances while the latter will in¬ 
duce strong ordering effects: Around a given ion, the 
first solvation shell will always be entirely constituted of 
oppositely charged species; this ordering automatically 
transfers up to several solvation shells. The third term 
of the intermolecular interactions, the dispersion, arises 
from correlated fluctuations of the electrons. It is always 
attractive and although it brings a smaller contribution 
than the other terms to the overall energy, it has impor¬ 
tant effects on the packing properties, thus influencing 
many thermodynamic properties such as the density or 
the surface tension. 

In a first approximation, the total energy arising from 
these three terms (repulsion, Coulomb and dispersion) 
can be written as a sum over all the pairs of ions. [HlfTS] 
From the simulation point of view, it is rather straightfor¬ 
ward to implement pair-potentials, so that the first sim¬ 
ulations of molten salts were performed relatively soon 
after the first molecular dynamics simulations were re¬ 
ported. [HjlfHil Significant results were obtained from 
these first simulations. For example the bulk proper¬ 
ties of several alkali halide crystals could be determined 
with a good precision, and a complete atomic-scale pic¬ 
ture of the structure and transport in simple molten salt 
was obtained for the first time. mm 

In a second stage, similar interaction potentials were 
tested for the study of more complicated halide melts in¬ 
cluding multiply charged cations, such as Zn 2+ in ZnCl 2 . 
The results were not so convincing; for example the pre- 
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dieted structures differed qualitatively from the exper¬ 
imental ones, which had been obtained from neutron 
or X-ray diffraction. [221425] Although the tetrahedral 
structure of the Zn 2+ first-neighbour coordination shell 
was well reproduced, the corner-sharing arrangement be¬ 
tween the tetrahedra was not predicted by the simula¬ 
tions, mum unless unphysically large dispersion inter¬ 
actions between the cations were introduced. [28] This 
inconsistency was solved by introducing a fourth compo¬ 
nent into the interaction potential, which represents the 
polarization of the electron cloud of an ion in response 
to its environment. This term cannot be imple¬ 

mented in the form of a pair-potential because of its in¬ 
herent many-body character. The simulations including 
polarization effects on the chloride anions were able to re¬ 
produce the structure factors extracted from experiments 
not only for pure ZnCl 2 but also for its mixtures with al¬ 
kali chlorides [3L, which allowed for an understanding 
of the basic mechanisms for network formation in ionic 
fluids. 

A systematic set of polarizable potentials for ionic 
solids based upon the Shell Model [33] had been proposed 
by Sangster and Dixon. [55] Yet these potentials suffered 
from numerical instability problems (the so-called polar¬ 
ization catastrophe) especially when applied to melts at 
high temperature and pressure conditions. Madden and 
Wilson proposed a different approach, [3f)! which also 
included an important short-range polarization effect. 
This term was introduced as a consequence of a study 
of polarization in condensed phases based upon first- 
principles electronic structure calculations [84] - rather 
than as a simple empirical stabilization device to avoid 
the polarization catastrophe. The model which was in¬ 
troduced in this study is now widely used in the study 
of ionic materials, ranging from molten fluorides |35l 136] 
to silicates [37] [35] and silver halides. [39, 40] By relat¬ 
ing the parameters of the interaction model to proper¬ 
ties of the individual ions (ion size, polarizability etc.) 
in a chemically-inspired way the model was shown able 
to account systematically for the solid state structures 
adopted by a very wide range of materials [44]. More 
recently, methods have been introduced to allow most of 
the parameters which enter into the interaction models 
to be derived directly from first-principles, rather than 
allowing them to be determined empirically by adjust¬ 
ing them to reproduce the experimental properties of 
the simulated material [42]. These ab zmfio-determined 
properties have included polarizabilities 45], the short- 
range repulsion and also the dispersion interactions 144] , 
so that the objective of finding realistic, predictive and 
transferable models has been accomplished, at least for 
halide melts. 

In the first section, we will introduce the correspond¬ 
ing functional form; with special emphasis given to the 
methodological aspects associated with the polarization 
effects. The second part details the consequences for the 


physico-chemical properties which relate to structure and 
dynamics. The interfacial properties, which are impor¬ 
tant in many industrial applications of ionic melts, are 
also discussed. We focus on the systems for which most 
of the data are available, i.e. inorganic ionic compounds, 
but the particular case of the room-temperature ionic 
liquids, a special class of molten salts in which at least 
one species is organic, will also be briefly discussed to 
indicate how the ideas gained from the study of “simple” 
molten salts are being transferred to these more complex 
systems. The ab initio parameterization of the potential 
parameters will not be considered in any detail here, it 
has been the subject of another recent review [42] . 


INTRODUCING POLARIZATION EFFECTS IN 
MOLECULAR DYNAMICS SIMULATIONS 

The focus of this topical review are the polarization 
effects themselves, so that other many-body effects that 
may occur in ionic liquids like the “breathing” of the 
anions will not be discussed. |45| Under this condition, 
the charge-charge, repulsion and dispersion contributions 
to the energy can be expressed through: 


/ n i n j .... (^ l 3 

^ RIM = E ( 7W + /sE ij ) 

(1) 

where the superscript RIM stands for “rigid ion model”. 

In the framework of an ionic model, which is based upon 
the interactions of closed-shell species, the charges q 1 
should be the formal, valence ones (except in the spe¬ 
cial case of redox active species, which has hardly been 
tackled by molecular dynamics simulations [46]). A*- 7 , 
a lJ , Cg and are parameters which have to be set up 
for each ion pair, and are Tang-Toennies dispersion 
damping functions, m describing the short-range pene¬ 
tration correction to the asymptotic multipole expansion 
of dispersion, [48] which take the following form: 
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Note that in many studies this correction is not taken into 
account, which is done by replacing the Tang-Toennies 
function by a constant value of 1.0 for any interatomic 
distance rW 

Compared to the RIM, the polarizable ion model 
(PIM) takes into account the effects arising from the po¬ 
larization of each ion: 
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This polarization, which is the response of the electron 
cloud of the ion i to the local electric field and its deriva¬ 
tives at the position r l , together with short-range effect 
to be described below, can be represented with a mul¬ 
tipole expansion. For example, up to the second order, 
there will be induced dipoles and quadrupoles on each 
ion: 

Ma as = ®apEp( r 1 ) + ^B a/ 3 75 ^(r 1 )F; 75 (r 1 ) + ... (4) 

= ^B a p 1 sE 1 {r 1 )Es(r') + C a p 1 sE 1 s{ r 1 ) + ... (5) 

where a and C are the dipole and quadrupole polarizabil¬ 
ities and B is the dipole-dipole-quadrupole hyperpolar¬ 
izability. IlHi E a and E a p are components of the electric 
field and field gradient, respectively. 

In the present manuscript, we will consider the case 
of dipole polarization effects only. As we shall see, in 
the case of many ionic species, this choice can safely be 
made without hindering the predictive capabilities of the 
model. In equation |4j the superscript as is used because 
it corresponds to the asymptotic contribution to the in¬ 
duced dipole only. Electronic structure calculations have 
showed that a substantial short-range effect on the in¬ 
duced dipole is not taken into account by this expres¬ 
sion. [34j !50j The origin of this short-range contribution 
to the dipole is schematized on figure [l] Compared to 
the perfect crystal cases where no induced dipoles are 
created, when a cation at a distance greater than next- 
nearest neighbour separation from the considered anion is 
displaced off its lattice site the latter feels an electric field 
which induces an asymptotic induced dipole /T ,as . But 
when a cation in the immediate vicinity of the anion is 
displaced off its lattice site, whilst the electric field tends 
to push the electrons away from the displaced cation 
(/T’ as ), the electrons also have more freedom to move into 
the space vacated (/T’ sr ). The total dipole induced by 
displacing first neighbour cations is reduced below that 
expected from the asymptotic term. On the contrary, 
when the central, polarizable species is the cation, the 
asymptotic and short-range dipoles generally point to¬ 
wards the same direction. |51 In practical applications, 
the short-range effects on the induced dipoles are taken 
into account by again introducing Tang-Toennies damp¬ 
ing functions m 
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where we draw the attention to the presence of an addi¬ 
tional parameter, cd, compared to equation [2] which is 
used for the dispersion interaction. This parameter mea¬ 
sures the strength of the ion response to the short-range 
effect and therefore depends on the identity of the ion. 
The efficiency of these functions in accounting for short- 
range effects was tested against ab initio calculations of 
the induced multipoles in distorted crystals. Enina 
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FIG. 1. Origin of the “asymptotic” and “short-range” con¬ 
tributions to the dipole induced on an anion in a crystal, a) 
The crystal is perfect, the anion does not feel electric field 
and there is no induced dipole, b) A cation at a distance 
greater than next-nearest neighbour separation from the con¬ 
sidered anion is displaced off its lattice site. The anion feels 
an electric field which induces an asymptotic induced dipole 
(i t ’ as . c) A cation in the immediate vicinity of the anion is 
displaced off its lattice site. Whilst the electric field tends 
to push the electrons away from the displaced cation (/T’“ s ), 
the electrons also have more freedom to move into the space 
vacated, which results in a short-range contribution to the 
induced dipole ( p, t,sr ). 
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The total induced dipole on each ion then satisfies the 
following equation: 

M* = a' E (<?iV iJ ')T (1 V - T (2 > ■ /**) (7) 

where we have introduced the charge-dipole and dipole- 
dipole interaction tensors: 

T«=vi =- \r ij (8) 

r ij 3 y ’ 

T (2) = V ® T (1) = ~^r ij 0 r ij - -Ei (9) 

r ij° r ij A 

and where I is the identity matrix. During a molecular 
dynamics simulation the induced dipoles have to be de¬ 
termined at each time step; the presence of the dipoles 
of all the ions other than i in the right hand term of 
equation [7] obliges us to solve the set of N equations 
self-consistently. As soon as the dipoles are known, the 
polarization energy is calculated from 

^polarization = E [(^VaSS^) “ VW 9d )) T a^°) 



where the first two terms arise respectively from the 
charge-dipole and the dipole-dipole interaction, and the 
last term corresponds to the the energy cost of deform¬ 
ing the charge density of the ion i (with polarizability a 1 ) 
to create the induced dipole. We immediately see that 
solving the set of equations [7] self-consistently is formally 
equivalent to solving: 

^polarization 
()fl' 

that is, to minimize ^polarization with respect to the in¬ 
duced dipoles. Note that, in the simulations, the po¬ 
larization effects which are intrinsically of a many-body 
nature, are calculated at the cost of evaluating only pair¬ 
wise additive interactions (plus the cost of enforcing the 
minimization condition)! The minimization task can 
be performed by using preconditioned conjugate gradi¬ 
ents methods, which appears to be much faster than the 
self-consistent approach. The dynamics is thus similar 
to the so-called Born-Oppenheimer ab initio molecular 
dynamics. As for the latter, algorithms based on the 
Car-Parrinello approach [55H55] or on predictor-corrector 
schemes HEn can also be used to determine the in¬ 
duced dipoles instead. 

The polarization potential (as well as the Coulombic 
one) consists of long-ranged interactions, for which it is 
necessary to go beyond the normal minimum image trun¬ 
cation and to sum the intermolecular interactions over 



the periodic images of the simulation cell. The most 
popular method for doing this is the Ewald summation 
technique, which consists in decomposing the interaction 
potential into a short-range component summed in real 
space and a long-range component summed in Fourier 
space. j53 (55] Several studies were devoted to the imple¬ 
mentation of this technique for multipole moments, 150 - 
l62j and the literature on this topic was recently reviewed 
by Stenhammar et al. [53] These authors pointed to sev¬ 
eral discrepancies between the formulae in the published 
work and they provided a consistent set of expressions. 

Expressions for the stress tensor and for other quan¬ 
tities like the heat current, which are obtained by dif¬ 
ferentiation of the expression for the total energy with 
respect to external variable, are also needed in order to 
perform molecular dynamics simulations in different en¬ 
sembles and to determine properties like the viscosity or 
thermal conductivity. For example, an expression for the 
stress tensor, required in NPT simulations or to calcu¬ 
late the viscosity, is obtained by differentiating the energy 
with respect to the shape of the simulation cell. [B3J [55] 
For a system in which the interaction between the par¬ 
ticles are described by short-range pair potentials, the 
stress tensor elements are given by: 

n^ = E m ’«-E^' (!2) 

i j>i ' p 


where m l is the mass of particle i. Here again, special 
care must be taken in the calculation of the long-ranged 
interactions, and the Ewald summation method must be 
employed. m In a simulation model which contains ad¬ 
ditional degrees of freedom, like the dipoles in the present 
polarizable potentials whose values will vary when the 
cell shape is changed, the question arises as to how they 
are to be handled in obtaining these expressions. The 
answer is that because the dipoles are in fact wholly de¬ 
termined by the instantaneous positions of all ions, be¬ 
cause of the minimization condition (equation 11), any 
additional terms due to the derivatives of dipoles with 
respect to cell-shape are to be ignored. The resulting 
expressions will still involve the dipole values and they 
must always be evaluated with the dipoles obeying the 
minimization condition. Once again, the expressions for 
the stress tensor etc. are then evaluated from only pair¬ 
wise additive expressions, despite the implicit many-body 
character of the interactions giving rise to them. 


ILLUSTRATION OF THE POLARIZATION 
EFFECTS 

Polarization effects on the structure 

The first ionic compounds which were systematically 
studied by molecular dynamics (and Monte-Carlo) simu- 
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FIG. 2. Partial radial distribution functions in pure LiF at 
1123 K. 


lations were the alkali halides. For such systems, the sim¬ 
ple RIM provides a rationalization of the crystal struc¬ 
tures depending on the two ionic radii. In all cases 
the unlike ion Coulombic interactions are maximized, 
and the maximal number of anions (cations) around a 
given cation (anion) is limited by their relative size only. 
Several typical structures can therefore be formed de¬ 
pending on their relative packing. For successively lower 
cation/anion radius ratio, these are the eight-coordinate 
caesium chloride (B2), six-coordinate rocksalt (Bl) and 
four-coordinated blende (B3) or wurtzite (B4) structures. 
The mechanisms involved in the successive transforma¬ 
tions when passing from one phase to another have also 
been investigated. [66] 

Upon melting, alkali halides keep a structure based 
on a competition between the overlap-repulsion and 
Coulombic interactions [13]. A typical set of radial dis¬ 
tribution functions (RDFs), obtained for pure LiF at 
1200 K, is provided on figure [2] The main feature to be 
noted is the position of the first maximum of the like-like 
RDFs, which occur at the same position coincident with 
the first minimum of the Li + -F _ one; the same behav¬ 
ior is also observed for all the following extrema. For all 
these systems, no structural change was observed when 
including the polarization effects. m 

As soon as ions with multiple charges are involved, the 
situation changes, there are multiple examples of multi¬ 
valent metal halide species for which polarization effects 
have shown to be important. Simulations based on the 
RIM model were unable to predict the structure of many 
MX 2 compounds (where M 2+ = Mg 2+ , Zn 2+ , Mn 2_t ", 
Ca 2+ , Sr 2+ or Ba 2+ and X - is an halide anion). This is 
particularly true for compounds with a low cation/anion 
radius ratio, which tend to crystallize in layered CdU or 
CdCU structures despite the fact that such structures 



FIG. 3. Snapshot of a representative configuration obtained 
for BeF 2 in the a-quartz structure. Blue: Be 2+ ions, Green: 
F _ ions, Red: Vector showing the direction and sense of the 
induced dipoles on F~ ions. Snapshot obtained using the 
VMD program. m 


involve shorter separations between the highly charged 
cations than could be achieved by other ways of accom¬ 
modating the cations in the anion lattice. 05] The role of 
polarization effects in the formation of these particular 
structures, which persist in the molten phase, has now 
been well identified. In fact all these structures depart 
from the picture provided by the RIM because they in¬ 
volve the formation of bent M-X-M angles, whereas the 
coulombic repulsion between the highly charged cations 
tends to push them as far apart as possible. The po¬ 
larization is the driving force for the occurrence of this 
non-trivial bond angle: When the halide anion is dis¬ 
placed off the line of centres of the cations a dipole is 
induced, which serves to screen the cation-cation repul¬ 
sion and lowers the total energy. Even for oxide materi¬ 
als, for which many-body effects other than polarization 
such as the aspherical breathing of the oxide anion O 2- 
have to be taken into account for an accurate and trans¬ 
ferable description [251 I58H70] . the polarization effects 
play a dominant role in determining these local ionic ar¬ 
rangements. 
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In order to illustrate this effect, we show in figure[3]the 
direction and sense of the induced dipoles on the F~ ions 
in the a-quartz structure of BeF 2 . It appears clearly that 
this dipole points towards the bisector of the Be-F-Be an¬ 
gle. In other words, the electronic cloud of the anion is 
shifted with respect to the position of the nuclei, which 
effectively screens the cation-cation Coulombic repulsion. 
In the framework of the RIM this screening could be mim¬ 
icked through the use of partial charges; this approach 
will be discussed in the next subsection, but we immedi¬ 
ately see that it suffers from transferability issues.When 
studying a family of ionic compounds such as molten flu¬ 
orides of varying compositions, the fluoride-fluoride in¬ 
teraction parameters should always be kept the same. 
This requirement holds for the Coulombic interaction pa¬ 
rameters, and this can only be achieved by using formal 
charges. Such a transferable rigid ion model (TRIM) was 
proposed by Woodcock et al. for a series of systems in¬ 
cluding KC1 and ionic liquids of MX 2 stoichiometry which 
have glass-forming ability (BeF 2 , ZnCl 2 , and Si0 2 ). [56] 
This model has proven very useful for understanding im¬ 
portant physical aspects of the pure materials, [MU] 
but the transferability of rigid ion models across a wide 
range of mixture compositions has never been proven. 

In recent years, a good deal of effort has been devoted 
to the study of molten fluorides because of their poten¬ 
tial use as a solvent in the molten salt fast reactor. im 
Experiments on molten fluorides are difficult because of 
the high melting points and their corrosive nature. A set 
of PIM interaction potentials has therefore been parame¬ 
terized on the basis of first-principles electronic structure 
calculations. m These potentials are therefore predictive 
in the sense that no empirical information is used in their 
construction. The parameterization is based on a gener¬ 
alized “force-matching” method. A suitable condensed- 
phase ionic configuration is taken from a molecular dy¬ 
namics simulation using some approximate force-field for 
the material of interest. Typically a hundred ions would 
be used in periodic boundary conditions. The configu¬ 
ration is then input to a planewave density functional 
theory (DFT) electronic structure program and an en¬ 
ergy minimization carried out to find the ground-state 
electronic structure. From the results of this calculation 
the force and dipole moment on each ion is obtained, the 
latter by making use of the transformation of the Kohn- 
Sham orbitals to a Maximally Localized Wannier Func¬ 
tion (MLWF) set. [771 The parameters in the polarizable 
potential are then optimized by matching the dipoles and 
forces from the potential on the same ionic configuration 
to the ab initio values. m If necessary the process may 
be iterated, by using the fitted potential to generate a 
new ionic configuration to input to the ab initio calcula¬ 
tion. In this approach, the dispersion terms have to be 
determined separately because of the use of functionals 
(e.g. PBE ,78]) for the DFT calculations in which these 
effects are not properly taken into account - Note that 


this difficulty will probably be overcome in the nearest 
future thanks to the development of improved function¬ 
als m The resulting potentials may be used in much 
larger scale molecular dynamics simulations to obtain the 
physical properties of interest. 

The procedure was validated on the LiF-BeF 2 mix¬ 
tures. This choice was guided by the existence of a 
comprehensive experimental database for these materi¬ 
als as a legacy of the Molten Salt Reactor Program con¬ 
ducted in the US in the 60’s. The predictive power of 
the simulations was tested against several experimental 
datasets without any empirical adjustments of the poten¬ 
tials. This included thermodynamic and transport data; 
at the atomistic structure level, the X-ray diffraction pat¬ 
terns m as well the infrared and Raman spectra [7611801 
were reproduced extremely well. The picture of a net¬ 
work of BeF 4 tetrahedral entities that are connected by 
their corner has emerged. Depending on the BeF 2 con¬ 
centration, the proportions of the various fluoroberyllate 
anions (BeF 2- , Be 2 Fi]! - , BeaF^g etc ) could be quanti¬ 
fied. I® [85] The important role of the induced dipoles 
in the transferability of the LiF-BeF 2 potential is high¬ 
lighted in figure [4] It shows the probability to find an 
induced dipole on a fluoride anion with respect to its 
intensity and to the angle that it forms with the corre¬ 
sponding Be-F bond. Surfaces appearing in a red color 
correspond to the highest probability whereas the white 
ones correspond to a probability of zero. Two sets of data 
have been separated: The top panel correspond to the 
bridging fluorides, which are shared by two Be 2+ atoms. 
This figure is in agreement with the picture provided on 
figure [3] (which is normal since all the fluoride anions are 
bridging ones in a-quartz BeF 2 ): The most likely angle 
is of around 120 degrees, i.e. when the induced dipole 
is directed along the Be-F-Be angle bisector. The bot¬ 
tom panel corresponds to the terminal fluoride, which are 
linked to one Be 2+ only. In that case the induced dipole is 
almost always directed along the Be-F bond, with a most 
likely angle of 160 degrees. The induced dipoles therefore 
allow the fluoride anions to adapt to distinct structural 
environments, which is not possible when using simpler 
pair-potentials only. 

The validation of the first-principles procedures on the 
LiF-BeF 2 mixtures has allowed us to adopt a predictive 
strategy for physico-chemical properties that hitherto re¬ 
main unknown despite their importance in establishing 
industrial proceses involving molten salts. This was done 
recently for ZrF 4 mixtures, for which a potential that was 
developed to predict thermodynamic properties [83] 184] 
was successfully used to interpret high-temperature EX- 
AFS spectroscopy data. [56] [SS] The structure of molten 
lanthanide chlorides, which involve trivalent cations, has 
also been elucidated thanks to the use of PIM potentials. 
Those were able to reproduce the neutron diffraction pat¬ 
terns [86H88] as well as Raman [TT] [85J and EXAFS [SU] 
spectrocopy data. In conclusion, it appears that taking 
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FIG. 4. Probability to find an induced dipole on a fluoride 
anion with respect to its intensity and to the angle that it 
forms with the corresponding Be-F bond. Surfaces appearing 
in a red color correspond to the highest probability whereas 
the uncoloured ones correspond to a probability of zero. Data 
obtained for molten Li 2 BeF 4 at a temperature of 873 K. 


the polarization effect into account is mandatory in any 
structural study of ionic materials involving multivalent 
species. 


On the use of partial charges in effective potentials 


As outlined above, in MX„ compounds the main role 
played by the polarization of the anion X - on the local 
structure is to screen the strong repulsive electrostatic 
term between two cation M” + in the coordination shell 
of the anion, which affects the M-X-M angle and may 
influence the connectivity of the cation-centred coordi¬ 
nation complexes (corner- or edge-sharing). Of course 
the drawback is an increase in the computational cost 
(mainly due to the self-consistent determination of the 



FIG. 5. Partial radial distribution functions in pure Z 11 CI 2 
at 1000 K. 


induced dipoles at each time step), so that it is tempting 
to find an ’’effective” way to account for it in the frame¬ 
work of the RIM. The simplest way to perform this is to 
reduce directly the cation-cation Coulombic interaction 
by using partial charges for the ions. This is routinely 
done in the case of molten silicates pH H2 and other 
oxide-based compounds. Most of the time, the param¬ 
eters of the potential (including partial charges) are fit¬ 
ted to reproduce bond length and angles, [93H95] so that 
these potentials can provide a good representation of the 
structure. We will show in the next section that it is 
much more difficult to simultaneously reproduce the dy¬ 
namical properties. The polarizabilities of the ions gives 
rise to a value for the high (optical) frequency dielectric 
constant (eoo) greater than unity and at large interionic 
separations. The effect is to reduce the strength of the ef¬ 
fective coulombic interactions by a factor of e^ 1 [SB]. It is 
noteworthy that the reduction in the charges required in 
the successful silicate pair potentials is much larger than 
would be anticipated from this effect alone (see below). 

Much less literature is available in the case of halide 
compounds. One of the best documented examples for 
which the polarization effects are important is molten 
ZnCl 2 - This system became the focus of attention when 
the first X-ray J23J and neutron m diffraction works 
appeared. Both studies showed the existence of inter¬ 
mediate range order in the liquid as manifested by the 
presence of a first sharp diffraction peak (FSDP) at a 
scattering vector k = 1 A -1 in the diffraction pattern. 
The first interpretation of their data by Biggin and En- 
derby led them to conclude that this FSDP was due to 
the Zn-Zn correlations. m Although this conclusion has 
been challenged, [55] it was confirmed recently by Zei- 
dler et al. who have remeasured carefully the full set of 
partial structure factors. [231 The main particularity of 













their extracted RDFs is the position of the first peak for 
the Zn-Zn partial function, which occurs at the same po¬ 
sition as the Cl-Cl one despite the larger charge on the 
cation. 

Here we will only concentrate on the molecular dy¬ 
namics simulation work which has been devoted to this 
compound. As already mentioned, the TRIM potential 
developed by Woodcock et al. was able to provide the 
correct Zn-Cl first-neighbour distances, together with a 
structure organized in a tridimensional network. When 
the neutron data became available, it was clear that the 
Zn-Zn distance was too long because of the lack of screen¬ 
ing of the Coulombic interactions. Gardner and Heyes 
tried to use several sets of partial charges without improv¬ 
ing the situation. m A RIM with an unphysically large 
dispersion attraction between the Zn cations does bring 
about the shift in the first peak of the Zn-Zn RDF, but 
at temperatures and densities where real ZnCl 2 is highly 
fluid, the simulated ”fluid” is glassy. Including the po¬ 
larization effects in a PIM with formal charges gives the 
correct RDFs for the melt, which are shown on figure [5] 
reproduces the intermediate-range order seen in the first 
sharp diffraction peak |55j, and also gives good dynam¬ 
ical properties. The agreement with experiment extends 
to mixtures with alkali chlorides [31] . 

Polarization effects on the dynamics 

In condensed phases, the dynamical processes with the 
fastest relaxation time are associated with the vibrations 
of the atoms. In crystals, vibrational properties may be 
studied from the phonon dispersion relationship while in 
amorphous compounds (glasses, liquids) infrared and Ra¬ 
man spectroscopies provide information on the character¬ 
istic vibrational motions. As we have already indicated, 
in many melts, especially those with polyvalent cations, 
relatively long-lived, quasi-molecular coordination com¬ 
plexes form and these give rise to well-resolved charac¬ 
teristic vibrational bands. EHain] On longer time scales, 
in the liquid state or in superionic solids, I100H102] the 
dynamics is characterized by the diffusion of ionic species 
in response to local chemical potential gradients. These 
transport properties are then quantified by a diffusion 
coefficient for each species present in the system. m 
The diffusion coefficients, conductivity and viscosity in 
the melt are strongly influenced by the lifetimes of the 
coordination complexes and by the degree to which these 
complexes are linked together to form a network. [52, 104] 
The inclusion of the polarization effects in the interac¬ 
tion potential is crucial for describing these complexes 
and their network-forming tendencies correctly. 

All the vibrational and diffusive properties can be 
rather straightforwardly extracted from a molecular dy¬ 
namics simulation. Even if the use of partial charges in 
rigid ion models may suffice to bypass the inclusion of an 


explicit polarization term in the study of structural prop¬ 
erties of ionic compounds, to also reproduce the dynamic 
properties with such a structurally optimized model is 
far more demanding. In order to illustrate this, we take 
the example of germania, GeC> 2 , which is a close struc¬ 
tural analog of silica. [ 10511106] For this system, Oeffner 
and Elliott have obtained a first set of parameters for a 
rigid ion model by using a two step procedure. [94] In a 
first stage, they calculated the ab initio potential energy 
surface of a Ge(OH) 4 cluster, on which they could fit one 
hundred possible set of potentials. Among those, they 
chose the potential that better reproduced the charac¬ 
teristic bond lengths and bond angles for the a-quartz 
GeC >2 structure. This potential involved for example 
partial charges of +1.5 e and -0.75 e for the Ge and 
O atoms respectively. It was used by Hawlitzky et al. 
1107] in simulations of the liquid state, and was shown to 
provide diffusion coefficients compatible with the exper¬ 
imental data (which consists in a series of viscosity mea¬ 
surements, from which the diffusion coefficients were es¬ 
timated using Eyring equation). Nevertheless, as soon as 
vibrational properties are concerned, important discrep¬ 
ancies have been observed. For example, the vibrational 
density of states predicted too high frequencies for all the 
bands found for a-quartz Ge02, which led Oeffner and 
Elliott to propose a second set of parameters, in which all 
the interactions were rescaled (for example, the new par¬ 
tial charge is of 0.94174 e for Ge atom). [M] This poten¬ 
tial now provided a better description of the vibrational 
properties, but then its use in liquid state simulations 
yielded values of the diffusion coefficients which were one 
order of magnitude too high. [10'0[ 1103 

Interaction potential parameters were obtained for a 
PIM (with formal valence charges) of Ge02 by Marroc- 
chelli et al. by force-fitting to first-principles DFT cal¬ 
culations m- The analytic expression for the repul¬ 
sion term used in that study differed slightly from the 
one given in equation [l] in order to enhance the stabil¬ 
ity of the potential against numerical problems in high 
temperature simulations. This potential was shown to 
provide diffusion coefficients in agreement with the orig¬ 
inal Oeffner-Elliott potential together with very good vi¬ 
brational properties. The latter could be tested by cal¬ 
culating the infrared spectrum in the glassy state. In 
fact, the inclusion of polarization effects for the oxide 
ions in the model may influence the predicted spectrum 
in two ways. um First, the interactions of the oxide 
ion dipoles may alter the local structure of the network 
and the strength of the bonds, which may introduce a 
shift of the vibrational frequencies. Second, the induced 
dipoles will themselves be responsible for absorption, as 
they too contribute to the total polarization fluctuations. 
The absorption coefficient in the presence of these extra 
moments is calculated from the imaginary part of the 
total dielectric function [n(^)a(i/) = 27w3 ! (e(i'))], which 
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can be determined following Caillol et al. mm mu as 

e{v) - £oo = ((M(0) 2 ) + 2niv(M • M)„ (13) 

3eo r 

+2(M-J) i , + ^-(J- J>„) 

where 

/»oo 

(3-3)„ = e 2mvt (3{t) ■ J(0))di, (14) 

Jo 


J(f) is the charge current J (t) = an d M(t) 

is the total system induced dipole moment, M(t) = 
JT/n*(t). The agreement was found to be very satis¬ 
factory, with a good correspondence of peak frequencies 
with experiment and a much better distribution of inten¬ 
sity across the spectrum than obtained by omitting the 


induced dipole terms in equation 14 



FIG. 6. Decomposition of the total infrared absorption spec¬ 
trum of glassy GeC >2 in its three component parts (JJ: charge- 
charge, MJ: charge-dipole, MM: dipole-dipole). 

The importance of polarization effects in determining 
the relative intensities of the bands can be demonstrated 
by separating the various contributions to the absorp¬ 
tion spectrum. In the case of glassy silica [Mj and 
beryllium fluoride EE it was observed that the interfer¬ 
ence between the induced dipoles and permanent charge 
contributions to the total polarization, contained in the 
(M • 3) v cross term, is responsible for the changes of the 
relative intensities. Figure [ 6 ] shows the decomposition of 
the absorption spectrum in the case of glassy GeC> 2 . The 
two main contributions are the charge fluctuation and 
the cross term; both show bands at the same intensities. 
Note that the cross-term strongly reduces the intensity 
of the two low frequency bands relative to that which 
would be obtained from the charge fluctuations alone, 
whereas for the high-frequency band the cancellation is 


much weaker. For the high frequency band the relation¬ 
ship of the charge-charge and cross-terms is different to 
that found previously for SiC >2 and BeF^ jTBl (109] . In the 
latter cases, the charge-charge and cross-terms had the 
same sign, so that the net band intensity was slightly 
larger than that predicted by the charge-charge term 
alone. This difference in the behaviour of the calculated 
spectra might arise from the inclusion of an anion-anion 
damping term in the polarization part of the interaction 
potential here, which was not the case for the previous 
studies. 

The ab initio -parameterized interaction potential was 
subsequently used to study the behavior of GeC >2 under 
pressure in the glassy and liquid states. [112l It could 
predict a smooth transition from a tetrahedral to a oc¬ 
tahedral network with a significant number of pentaco- 
ordinated germanium ions appearing over an extended 
pressure range, in agreement with the most recent exper¬ 
imental data. m 

Polarization effects also influence the long-time dy¬ 
namics. Even in alkali halides, where the structural con¬ 
sequences of including polarization effects are not large, 
diffusion coefficients tend to be increased. This has been 
shown, for example, in studies of mixtures of LiCl and 
KC1 1671 of various compositions which were simulated by 
using both a rigid and a polarizable ion model in which 
all the repulsion, dispersion and ionic charges were kept 
exactly the same and, in the PIM case, a realistic polar¬ 
izability of the Cl~ ion was used. The comparison of the 
two sets of diffusion coefficients showed that the effect 
of including polarizability is to make the system more 
mobile. Interestingly, the largest increase concerns the 
Li + ions, although these ions are not polarizable. The 
calculated diffusivities for the pure melts appeared to be 
much closer to the experimental values for the PIM. Sim¬ 
ilar calculations were performed for a room temperature 
ionic liquid, namely the l-ethyl-3-methylimidazolium ni¬ 
trate. 116 . j ll7] Here again, the effect of including polar¬ 
ization effects was to enhance the fluidity of the liquid. 
This reflected on the whole set of transport coefficients: 
With the polarizable model the diffusion coefficients and 
ionic conductivities were increased, and consistently the 
viscosity was decreased. Far more dramatic increases in 
fluidity are seen in systems with polyvalent cations, here 
the issue is complicated by making a direct comparison 
between a PIM and an appropriate RIM. In many cases, 
simply omitting the polarization terms from a PIM gives 
an RIM which has hopeless structural and thermody¬ 
namic properties (which is not the case for alkali halides). 
It seems better to compare a PIM with an independent 
RIM which has been optimized, so far as possible, to re¬ 
produce the structure. In simulations of UCI 3 , m for 
example, such a comparison showed order of magnitude 
increases in the diffusivities of both ions with the PIM re¬ 
producing the measured conductivity and viscosity very 
well. 
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System T (K) A RIM A PIM A exp 
LiCl 1200 0.841 0.643 0.534 

NaCl 1300 0.581 0.509 0.478 

KC1 1300 0.387 0.343 0.345 

TABLE I. Values of the thermal conductivity (in Wm _1 K _1 ) 
for a series of molten chlorides obtained using PIM and RIM 
interaction potentials, and from experiments. m 

From the mechanistic point of view, these results are 
well understood in terms of the screening of the charge- 
charge interactions by the induced dipoles. The polariza¬ 
tion energy always lowers the total coulombic energy - it 
is calculated by minimizing the latter (equation [ll]). The 
effect is always largest for ionic configurations in which 
the ions are most strongly polarized, which is the case 
when they sit in an asymmetric coordination environ¬ 
ment. If one visualizes diffusion as ions hopping between 
locally metastable configurations, it is easy to imagine 
that the transition states in such hops are more asym¬ 
metrically coordinated and hence that polarization will 
lower the energy more at the barrier tops than at the 
minima. 

The picture of polarization effects lowering a barrier to 
structural relaxation accounts for the influence on those 
transport coefficients, like diffusion or viscosity, which 
depend on such events. The situation is somewhat more 
complicated for the thermal conductivity which measures 
a material’s ability to conduct heat, and is given, for a 
binary mixture, by m 

x=t ~ 2 { Lee ^t^)- (15) 

where 

Lab = Wk~B h ( 16 ) 

For systems consisting of more than three charged 
species, different expressions need to be derived for 
A. [120] Here the quantities involved are the charge cur¬ 
rent (a = Z), which is defined as previously, and the 
energy current (a = E). The computation of the lat¬ 
ter quantity requires to take special care due to the use 
of the Ewald summation technique. The first expression 
was derived by Bernu and Vieillefosse in their study of the 
transport coefficients of the one-component plasma, m 
and we extended this work to the case of potentials in¬ 
cluding polarization effects. m 

The values obtained for a series of molten chlorides 
(LiCl, NaCl, KC1) in the latter study 122; are sum¬ 
marized in table [I] In this work, the PIM parameters 
were obtained wholly from first-principles DFT calcula¬ 
tions, and the RIM simply corresponds to the same po¬ 
tential where the polarization effects are omitted. We 


observe that the PIM systematically predicts lower ther¬ 
mal conductivities, which are in much better agreement 
with the experimental data. [123] The non-polarizable 
model yields similar values to the famous Tosi-Fumi po¬ 
tentials, IHUI5] HU] which are known to give the correct 
structural properties of the MX molten salts in general. 
This shows again that the inclusion of polarizability is 
crucial for a good representation of the thermal relax¬ 
ation. The thermal conductivity is affected by different 
aspects of the ionic dynamics than the other transport 
coefficients, such as the viscosity or ionic conductivity; 
this quantity can therefore be used to test the interac¬ 
tion model in different ways. For the set of alkali halides, 
the first-principles-deterniined PIM potentials predicted 
values for all transport properties with the largest dis¬ 
crepancy being about 10% of the measured value. [1221 

Polarization effects on the interfacial properties 

Liquid-vapor interface 

The simplest interface involving an ionic liquid that 
one may consider is the one with its own vapor. The ex¬ 
perimental vapor pressure of these systems is very low, 
which means that from the point of view of the molec¬ 
ular simulation this interface is in reality an ionic liquid 
- vacuum one. The simulation cell now has dimension 
L x L x L z , where L z = D + L V acuum- Due to the use of 
vacuum periodic boundary conditions in the Ewald sum¬ 
mation instead of conducting boundary conditions for 
the calculation of long-range electrostatic interactions, 
additional terms have to be added to the energy, forces 
and pressure tensor. m It is also necessary to perform 
an Ewald summation of dispersion interactions, [126] to 
avoid substantial truncation effects. m 

The electrostatic potential difference across the inter¬ 
face contains two contributions, 

A$(s) = A$,(z) + AQ^z) (17) 

= -(-/' d/ f p q (z")dz" + f Puiz'jcUm 
eo y Jzo Jzo JZ 0 J 

which respectively correspond to the distribution of 
charges and induced dipoles. In this expression Zo corre¬ 
sponds to a point in the vapor region. The importance of 
the two terms is showed on figure [7] which corresponds to 
molten LiBeF 3 at a temperature of 1060 K. [128] In this 
system, we observed the stabilization of the fluoroberyl- 
late species at the interface, leading to a local enhance¬ 
ment of the Be 2+ ions concentration. This segregation 
leads to important charge-separation effects at the inter¬ 
face, which explains the very large value of the charge 
term. An electronic dipole moment is also created on the 
fluoride atoms at the interface, which results in a contri¬ 
bution opposite in sign, thus softening the electrostatic 
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FIG. 7. Electrostatic potential profile across the liquid-vapor 
interface in molten LiBeFa at 1060 K. The contributions due 
to the charges and to the induced dipoles are also shown. 


potential variation. The total electrostatic potential dif¬ 
ference between the vacuum and the ionic liquid is of 
approximately 0.98 V. 

Although it appears clearly that without the inclusion 
of polarization effects the structure of the liquid - vacuum 
interface of LiBeFa would differ a lot, no explicit compar¬ 
ison has been made for this system. Such a comparison of 
the RIM and PIM interaction potentials was performed 
in the case of KI, in a series of papers dedicated to a bet¬ 
ter understanding of the interfacial properties of molten 
salts. [12711129L1130J It was shown that polarization tends 
to increase the interfacial width by reducing slightly the 
bulk density and enhancing the ion density in the interfa¬ 
cial region. Stronger effects were observed for the surface 
tension, which can be obtained through the mechanical 
definition 

7 = ^(n 2Z ) — -(n xx + n yy )^ (19) 

where the various Il aa are the diagonal pressure tensor 
components. In the case of KI, the effect of polariza¬ 
tion is to reduce the surface tension by about 20 %. By 
calculating the Coulombic, repulsion, dispersion and po¬ 
larization components for the surface tension, Aguado et 
al. could show that the polarization component is small 
and positive, so that the effect of polarization on 7 is indi¬ 
rect: It is a consequence of the increase of the interfacial 
width. |127| It is worth mentioning that the importance 
of polarization effects in the behavior of ions at the air 
- water interface has also been the object of numerous 
studies. [1311 - 11 "33] 


The ionic liquid - electrified metal interface 

The study of the interface which is formed when an 
ionic liquid is in contact with a metal is of particular im¬ 
portance; the interest derives from several sources. In the 
case of high temperature molten salts, understanding the 
corrosion mechanisms is essential for the design of Gen¬ 
eration IV nuclear reactors m\mm as well as for 
the improvement of pyrometallurgical applications, such 
as the Hall-Heroult process - which is the major indus¬ 
trial process for the production of aluminium. 11371 1138] 
As for room-temperature ionic liquids, they are used as 
electrolytes in electrochemical supercapacitors, [139H141] 
batteries, and fuel cells. ]l42j In all these applications, a 
better understanding of the properties of the ionic liquid 
in the vicinity of the charged surface would be beneficial. 

The study of such interfaces by molecular dynamics 
simulations has been made possible with the introduc¬ 
tion of a method to represent the polarization of a model 
metallic electrode which is maintained at a controlled 
electric potential difference. The system is period¬ 
ically replicated in the plane parallel to the electrodes. 
The metallic, constant potential condition is attained by 
minimizing a suitable energy function with respect to 
variable charges on the electrode atoms, following a pro¬ 
cedure suggested by Siepmann and Sprik in a different 
context. 11441 

A suitable interaction potential was obtained, on the 
basis of small scale DFT calculations, for a system con¬ 
sisting of a molten LiCl electrolyte and a solid aluminium 
electrode. In addition to the four terms already present 
in the PIM, it was found necessary to include an addi¬ 
tional short-range term between the metal atoms and the 
ions, 

N M 

V SR = ^ R e~ a ^ R{rij ~ d)2 r ij ■ li\ (20) 

i=i j= 1 

where r y is a unit vector, in order to obtain a satisfactory 
representation of DFT-calculated dipoles and forces (N 
and M respectively are the number of ions in the liquid 
and of atoms in the metal) for the ions close to the elec¬ 
trode surface. [45, This term is directed along the inter¬ 
atomic separation; it behaves like an electric field which 
distorts the charge cloud of a melt ion due to overlap- 
mediated interactions with the electrode atoms. The 
functional form in equation [ 20 ] is somewhat arbitrary, 
as the range of atom-ion separations at which these in¬ 
teractions are sampled, in the immediate vicinity of the 
electrode surface, is small and insufficient to determine 
the shape of the full function. This form allows the forces 
and dipoles to be adequately represented in the physically 
significant regions without bad computational problems 
in regions which are physically irrelevant because of the 
repulsive interactions between electrode and ions. 
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FIG. 8. Variation of the metal surface charge with the po¬ 
tential drop across the electrified aluminium - molten LiCl 
interface. 


The model was then used in simulations of this sys¬ 
tem for various values of the electrical potential differ¬ 
ence applied between the electrodes. The first relevant 
quantity that can be extracted from such simulations is 
the electrostatic potential profile. It differs from the one 
shown in figure [7] because the potential in the electrodes 
is not 0 V anymore, unlike for the liquid-vapor interface. 
It is also possible to calculate the average charge of the 
electrodes a. Then the potential of zero charge (PZC) 
corresponds to the potential drop across the electrode - 
liquid interface for which the accumulated charge on the 
electrode surface is zero on average. 

A<f>pzc = A<f>0 = 0) (21) 


The calculated PZC cannot be compared to experimen¬ 
tally measured quantities since the latter are always mea¬ 
sured with respect to a reference electrode whereas we are 
dealing with absolute potentials. The variation of the 
surface charge with respect to Ad> was calculated in the 
case of the electrified aluminium - molten LiCl interface 
with and without the inclusion of polarization effects for 
the chloride anions (the lithium cations are not polar¬ 
izable). (1451 Three important differences are observed 
when polarization effects are included: First, the PZC 
goes toward a more negative value. Second, the differen¬ 
tial capacitance, which is given by the derivative of this 
function, 


C = 


da 

(9A*I> 


( 22 ) 


becomes higher, switching from 6.1 /rF cm -2 to 
10.2 /iF cm -2 . The analysis of the induced dipoles dis¬ 
tribution showed that they bring an additional screening 
of the electrode charge, 1145 : and that the magnitude of 
the dipole increases with the electrode potential. 


But the most important difference between the two 
sets of data is the presence of a discontinuity when the 
polarization effects are included. This increase is due 
to a potential-driven ordering transition in the adsorbed 
layer of electrolyte, which occurs at the PZC. The ad¬ 
sorbed layer switches from a fluid-like structure to a 2- 
dimensional crystalline one which is formed through an 
epitaxial mechanism to adapt to the electrode surface 
structure. This crystalline structure did not form when 
the RIM was used, which shows again that taking polar¬ 
ization effects into account is obligatory when studying 
the interfacial properties of molten salts - even for an 
alkali halide, where the polarization effects on the struc¬ 
ture of the melt itself are very small. For RTILs, the 
situation is much more complex: IT46I 1147 ) The internal 
charge distribution within the molecular ions already in¬ 
duces additional screening of the electrode charges and 
electronic polarization effects may play a secondary role. 
For example, a simple RIM was able to yield an or¬ 
dered surface structure for the adsorption of l-butyl-3- 
methylimidazolium hexafluorophosphate on a graphite 
electrode. m 


CONCLUSIONS AND PERSPECTIVES 

In this topical review we have given an overview of 
the importance of polarization effects on the physico¬ 
chemical properties of ionic materials, with particular 
focus on their melts. Most of the information was ex¬ 
tracted from molecular simulation studies, in which the 
choice of the model ( e.g. the rigid ion model versus the 
polarizable ion model) will determine whether these ef¬ 
fects are taken into account. In most halide and oxide 
materials, the most polarizable species is the anion. The 
creation of induced dipoles affects the thermodynamics of 
the system in different ways. From the structural point of 
view, although polarization does not seem to play an im¬ 
portant role for simple alkali halides, it becomes very im¬ 
portant as soon as multiply charged cations are involved. 
There, the induced dipoles screens the repulsive cation- 
cation Coulombic interaction, thus allowing the forma¬ 
tion of bent M-X-M angles and thereby influencing the 
character of the networks which form through the linkage 
of cation-centred coordination complexes . Most of the 
attempts to introduce this screening effectively by using 
partial charges failed in predicting the correct structural 
properties. Concerning the dynamic properties, polariza¬ 
tion seems to systematically increase the fluidity of the 
liquids, which results in higher diffusion coefficients for 
all the ions (including the non-polarizable ones such as 
Li + ) and the electrical conductivity as well as in a lower 
viscosity. Finally, the recent work devoted to the study of 
the interfacial properties of molten salts have shown the 
importance of the induced dipoles for the stabilization 
of particular surface structures such as the formation of 
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an ordered layer of LiCl on the (100) surface of metallic 
aluminium. 

Alongside the various qualitative points discussed 
within this manuscript we have also emphasized the im¬ 
portance of reliable interaction potentials. A substantial 
effort was devoted to developing such potentials in re¬ 
cent years. [42] The polarizable ion model, because it 
introduces additional degrees of freedom which mimic 
the response of the electronic structure of the ions to 
their changing coordination environments, can routinely 
be parameterized by fitting the predicted forces and 
dipoles to a large body of information generated from 
first-principles calculations. The resulting potentials are 
of first-principles accuracy, which means that molecular 
dynamics of ionic compounds can now be used as a pre¬ 
dictive tool under various conditions. 

Yet it is important to notice that there is still room 
for improvement of the models. For example, some en¬ 
vironmental effects occur on the electronic structure of 
the ions when the nature of the liquid is changed. In flu¬ 
oride molten salts, for example, calculations of the con¬ 
densed phase ionic polarizabilities have shown that the 
fluoride anion polarizability may shift a lot, passing from 
7.8 a.u. in molten LiF to 11.8 a.u. in CsF. }l49l This is 
due to differences in the confining potential, which affects 
the electron density around a given anion, and originates 
from both Coulombic interactions and the exclusion of 
electrons from the region occupied by the electron den¬ 
sity of the first-neighbor shell of cations. [M, 1501 151 j 
When passing from one cation to another (for example 
in the series Li + —>Na + —>K + —tCs + ), two effects are then 
competing: On the one hand, the anion-cation distance 
increases, which results in a diminution of the confining 
potential, but on the other hand, the volume occupied by 
the cation electron density also increases, with an oppo¬ 
site effect on the confining potential. Here, the observed 
increase of polarizability with the size of the cation tends 
to show that the first effect is the most important. In¬ 
deed, the value obtained for CsF is approaching the free 
F~ anion polarizability, which is 16 a.u. These effects 
have also been observed in similar calculations performed 
on solid oxides 03 an d protic solvents [15211153] . Such 
a result means that in order to build a completely trans¬ 
ferable interaction potential, the model should allow the 
variation of the polarizability in response to the fluctu¬ 
ations of the environments. First steps have been made 
in that direction in the case of MgO only. (BSJ EE. 

In the future it is probable that the polarizable 
force fields will also systematically be used for room- 
temperature ionic liquids. Although several studies have 
recently been reported, [TTB], 11T7 ; I154H158] the general 
approach of most groups remains to discard polarization 
effects. In fact, the situation seems to be rather similar 
to the case of simple alkali halides, with a correct rep¬ 
resentation of the structural properties and a poor esti¬ 
mation of the transport properties in general. nm urn 



FIG. 9. Variation of the electronic density when an exter¬ 
nal field is applied to an isolated l-ethyl-3-methylimidazolium 
cation. Green and red zones respectively correspond to an in¬ 
crease (decrease) of the electronic density with respect to the 
situation without the external field. 


From the technical point of view, the main difference 
between inorganic salts and room temperature ionic liq¬ 
uids arises from the molecular structure of the latter. As 
depicted on figure [9] where we plot the variation of the 
electronic density when an external field is applied to an 
isolated l-ethyl-3-methylimidazolium cation, there now 
are some intramolecular polarization effects that have to 
be taken into account. This can be done efficiently using 
the approach proposed by Thole, jl'S'D] which has now 
been introduced in numerous molecular dynamics simu¬ 
lation codes. 
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